﻿using OSGeo.GDAL;
using OSGeo.OGR;
using OSGeo.OSR;
using System;
using System.Collections.Generic;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
using System.Linq;
using System.Runtime.InteropServices;

namespace ConsoleApp
{
    public class MultiGeoTiffSplitter
    {
        public const string GEO_TIFF_SPLIT_DIR_NAME = "tiffsplits";
        public const string IR_GEO_TIFF_SPLIT_MAP_JSON_FILE_NAME = GEO_TIFF_SPLIT_DIR_NAME + "_T.json";
        public const string VS_GEO_TIFF_SPLIT_MAP_JSON_FILE_NAME = GEO_TIFF_SPLIT_DIR_NAME + "_V.json";

        private static GeoPoint PixelToGeo(int x, int y, double[] transform,
                                    CoordinateTransformation ct)
        {
            double geoX = transform[0] + x * transform[1] + y * transform[2];
            double geoY = transform[3] + x * transform[4] + y * transform[5];
            double[] argout = new double[6];
            ct.TransformPoint(argout, geoX, geoY, 0);
            return new GeoPoint { Lon = argout[1], Lat = argout[0] };
        }

        public static void SplitMultipleGeoTiffs(List<string> inputPaths, List<string> types,
                                     string outputDir, int tileSize = 512, bool overwrite = false)
        {
            if (inputPaths.Count == 0 || inputPaths.Count != types.Count)
                throw new ArgumentException("输入文件路径和类型数量不匹配");

            Directory.CreateDirectory(outputDir);
            GdalConfiguration.ConfigureGdal();
            GdalConfiguration.ConfigureOgr();
            Gdal.AllRegister();
            Ogr.RegisterAll();

            // 存储所有文件的地理信息和数据集
            List<FileGeoInfo> fileInfos = new List<FileGeoInfo>();
            List<Dataset> datasets = new List<Dataset>();

            try
            {
                // 1. 收集所有文件的地理信息
                foreach (string inputPath in inputPaths)
                {
                    if (!File.Exists(inputPath))
                        throw new FileNotFoundException("输入文件不存在", inputPath);

                    Dataset ds = Gdal.Open(inputPath, Access.GA_ReadOnly);
                    if (ds == null)
                        throw new InvalidOperationException($"无法读取GeoTIFF文件: {inputPath}");

                    datasets.Add(ds);
                    double[] geoTransform = new double[6];
                    ds.GetGeoTransform(geoTransform);

                    // 验证地理变换参数
                    if (Math.Abs(geoTransform[2]) > 0.001 || Math.Abs(geoTransform[4]) > 0.001)
                        throw new NotSupportedException("当前不支持旋转坐标系");

                    SpatialReference srcSrs = new SpatialReference(ds.GetProjection());
                    SpatialReference destSrs = new SpatialReference("");
                    destSrs.SetWellKnownGeogCS("WGS84");
                    CoordinateTransformation ct = new CoordinateTransformation(srcSrs, destSrs);

                    // 计算实际地理范围（考虑Y方向翻转）
                    double minX = geoTransform[0];
                    double maxX = geoTransform[0] + ds.RasterXSize * geoTransform[1];
                    double maxY = geoTransform[3];
                    double minY = geoTransform[3] + ds.RasterYSize * geoTransform[5];

                    FileGeoInfo info = new FileGeoInfo
                    {
                        Dataset = ds,
                        Transform = geoTransform,
                        Ct = ct,
                        XSize = ds.RasterXSize,
                        YSize = ds.RasterYSize,
                        MinX = minX,
                        MaxX = maxX,
                        MinY = minY,
                        MaxY = maxY,
                        PixelSizeX = geoTransform[1],
                        PixelSizeY = Math.Abs(geoTransform[5])  // Y方向像素大小取正值
                    };

                    fileInfos.Add(info);
                }

                // 2. 计算所有文件的公共重叠区域
                double commonMinX = fileInfos.Max(f => f.MinX);
                double commonMaxX = fileInfos.Min(f => f.MaxX);
                double commonMinY = fileInfos.Max(f => f.MinY);
                double commonMaxY = fileInfos.Min(f => f.MaxY);

                if (commonMinX >= commonMaxX || commonMinY >= commonMaxY)
                    throw new InvalidOperationException("文件之间没有重叠区域");

                // 3. 计算公共网格参数（使用所有文件的平均分辨率）
                double avgPixelSizeX = fileInfos.Average(f => f.PixelSizeX);
                double avgPixelSizeY = fileInfos.Average(f => f.PixelSizeY);

                // 4. 创建网格坐标系
                GridCoordinateSystem gridSystem = new GridCoordinateSystem(
                    commonMinX, commonMaxY, // 左上角
                    avgPixelSizeX, avgPixelSizeY
                );

                // 5. 计算网格尺寸
                int gridWidth = (int)Math.Round((commonMaxX - commonMinX) / avgPixelSizeX);
                int gridHeight = (int)Math.Round((commonMaxY - commonMinY) / avgPixelSizeY);

                // 6. 处理每个文件
                for (int fileIdx = 0; fileIdx < inputPaths.Count; fileIdx++)
                {
                    string type = types[fileIdx];
                    FileGeoInfo fileInfo = fileInfos[fileIdx];
                    List<GeoTileInfo> tiles = new List<GeoTileInfo>();

                    // 获取当前文件的颜色表
                    ColorTable colorTable = null;
                    using (Band firstBand = fileInfo.Dataset.GetRasterBand(1))
                    {
                        colorTable = firstBand.GetRasterColorTable()?.Clone();
                    }

                    // 7. 在重叠区域内分块处理
                    for (int gridY = 0; gridY * tileSize < gridHeight; gridY++)
                    {
                        for (int gridX = 0; gridX * tileSize < gridWidth; gridX++)
                        {
                            // 计算当前网格的地理范围
                            double tileMinX = commonMinX + gridX * tileSize * avgPixelSizeX;
                            double tileMaxY = commonMaxY - gridY * tileSize * avgPixelSizeY;
                            double tileMaxX = tileMinX + tileSize * avgPixelSizeX;
                            double tileMinY = tileMaxY - tileSize * avgPixelSizeY;

                            // 转换为当前文件的像素坐标
                            var pixelCoords = ConvertGeoToPixel(
                                tileMinX, tileMaxY,
                                fileInfo.Transform,
                                fileInfo.Dataset.RasterXSize,
                                fileInfo.Dataset.RasterYSize
                            );

                            int xStart = pixelCoords.Item1;
                            int yStart = pixelCoords.Item2;

                            // 确保在图像范围内
                            if (xStart < 0 || yStart < 0 ||
                                xStart >= fileInfo.Dataset.RasterXSize ||
                                yStart >= fileInfo.Dataset.RasterYSize)
                            {
                                continue;
                            }

                            ProcessOverlapTile(
                                fileInfo.Dataset,
                                xStart, yStart,
                                tileSize,
                                outputDir,
                                fileInfo.Transform,
                                fileInfo.Ct,
                                colorTable,
                                overwrite,
                                type,
                                tiles,
                                gridX,
                                gridY,
                                tileMinX,
                                tileMaxY,
                                tileMaxX,
                                tileMinY
                            );
                        }
                    }

                    // 8. 保存元数据
                    File.WriteAllText(
                        Path.Combine(outputDir, type == "T" ? IR_GEO_TIFF_SPLIT_MAP_JSON_FILE_NAME : VS_GEO_TIFF_SPLIT_MAP_JSON_FILE_NAME),
                        Newtonsoft.Json.JsonConvert.SerializeObject(tiles, Newtonsoft.Json.Formatting.Indented)
                    );
                }
            }
            finally
            {
                // 确保释放所有数据集
                foreach (var ds in datasets)
                {
                    ds.Dispose();
                }
            }
        }

        // 将地理坐标转换为像素坐标
        private static Tuple<int, int> ConvertGeoToPixel(double geoX, double geoY,
                                                      double[] transform,
                                                      int rasterXSize, int rasterYSize)
        {
            // 使用逆地理变换公式
            double det = transform[1] * transform[5] - transform[2] * transform[4];

            if (Math.Abs(det) < 1e-10)
                throw new InvalidOperationException("地理变换矩阵不可逆");

            double invDet = 1.0 / det;

            // 计算逆矩阵
            double inv1 = transform[5] * invDet;
            double inv2 = -transform[2] * invDet;
            double inv4 = -transform[4] * invDet;
            double inv5 = transform[1] * invDet;

            // 应用逆变换
            double x = (geoX - transform[0]) * inv1 + (geoY - transform[3]) * inv2;
            double y = (geoX - transform[0]) * inv4 + (geoY - transform[3]) * inv5;

            // 返回整数像素坐标（四舍五入）
            return Tuple.Create(
                (int)Math.Round(Math.Max(0, Math.Min(rasterXSize - 1, x))),
                (int)Math.Round(Math.Max(0, Math.Min(rasterYSize - 1, y)))
            );
        }

        private static void ProcessOverlapTile(
            Dataset dataset,
            int xStart, int yStart, int tileSize,
            string outputDir,
            double[] geoTransform, CoordinateTransformation ct,
            ColorTable colorTable, bool overwrite, string type,
            List<GeoTileInfo> tiles, int gridX, int gridY,
            double tileMinX, double tileMaxY,
            double tileMaxX, double tileMinY)
        {
            int width = Math.Min(tileSize, dataset.RasterXSize - xStart);
            int height = Math.Min(tileSize, dataset.RasterYSize - yStart);
            int bandCount = dataset.RasterCount;

            // 跳过不完整的边缘块
            if (width < tileSize || height < tileSize)
                return;

            // 统一命名规则：基于网格位置和文件类型
            string fileName = $"geotiff_{gridX}_{gridY}_{type}.png";
            string filePath = Path.Combine(outputDir, fileName);

            if (!overwrite && File.Exists(filePath))
            {
                Console.WriteLine($"跳过 {fileName}");
                return;
            }

            // 读取所有波段数据
            byte[][] tileData = new byte[bandCount][];
            for (int b = 0; b < bandCount; b++)
            {
                using (Band band = dataset.GetRasterBand(b + 1))
                {
                    tileData[b] = ReadBandData(band, xStart, yStart, width, height);
                }
            }

            // 检查是否为空白切片
            if (IsTileBlank(tileData))
            {
                Console.WriteLine($"跳过空白切片 {fileName}");
                return;
            }

            // 颜色处理决策
            if (bandCount == 1 && colorTable != null)
            {
                SaveWithColorTable(tileData[0], colorTable, width, height, filePath);
            }
            else if (bandCount >= 3)
            {
                SaveAsTrueColor(tileData[0], tileData[1], tileData[2], width, height, filePath);
            }
            else
            {
                GeneratePseudocolor(tileData[0], width, height, filePath);
            }

            // 记录元数据
            tiles.Add(new GeoTileInfo
            {
                FileName = fileName,
                X = gridX,
                Y = gridY,
                TypeId = type,
                TopLeft = new GeoPoint { Lon = tileMinX, Lat = tileMaxY },
                BottomRight = new GeoPoint { Lon = tileMaxX, Lat = tileMinY },
                Center = new GeoPoint
                {
                    Lon = (tileMinX + tileMaxX) / 2,
                    Lat = (tileMinY + tileMaxY) / 2
                }
            });
        }

        private static byte[] ReadBandData(Band band, int xStart, int yStart,
                                         int width, int height)
        {
            byte[] buffer = new byte[width * height];
            band.ReadRaster(xStart, yStart, width, height,
                          buffer, width, height, 0, 0);
            return buffer;
        }

        private static void SaveWithColorTable(byte[] data, ColorTable ct,
                                             int width, int height, string path)
        {
            using (Bitmap bmp = new Bitmap(width, height, PixelFormat.Format8bppIndexed))
            {
                // 设置调色板
                ColorPalette palette = bmp.Palette;
                for (int i = 0; i < ct.GetCount(); i++)
                {
                    ColorEntry ce = ct.GetColorEntry(i);
                    palette.Entries[i] = Color.FromArgb(ce.c4, ce.c1, ce.c2, ce.c3);
                }
                bmp.Palette = palette;

                // 快速写入数据
                BitmapData bmpData = bmp.LockBits(new Rectangle(0, 0, width, height),
                                              ImageLockMode.WriteOnly,
                                              PixelFormat.Format8bppIndexed);
                Marshal.Copy(data, 0, bmpData.Scan0, data.Length);
                bmp.UnlockBits(bmpData);
                bmp.Save(path, ImageFormat.Png);
            }
        }

        private static void SaveAsTrueColor(byte[] r, byte[] g, byte[] b,
                                          int width, int height, string path)
        {
            using (Bitmap bmp = new Bitmap(width, height, PixelFormat.Format24bppRgb))
            {
                BitmapData bmpData = bmp.LockBits(new Rectangle(0, 0, width, height),
                                                ImageLockMode.WriteOnly,
                                                PixelFormat.Format24bppRgb);

                unsafe
                {
                    byte* ptr = (byte*)bmpData.Scan0;
                    for (int i = 0; i < width * height; i++)
                    {
                        ptr[2] = r[i]; // R
                        ptr[1] = g[i]; // G
                        ptr[0] = b[i]; // B
                        ptr += 3;
                    }
                }
                bmp.UnlockBits(bmpData);
                bmp.Save(path, ImageFormat.Png);
            }
        }

        public static void GeneratePseudocolor(byte[] data, int width,
                                              int height, string path)
        {
            // 自动计算极值
            byte min = byte.MaxValue, max = byte.MinValue;
            foreach (byte val in data)
            {
                if (val < min) min = val;
                if (val > max) max = val;
            }

            using (Bitmap bmp = new Bitmap(width, height, PixelFormat.Format24bppRgb))
            {
                BitmapData bmpData = bmp.LockBits(new Rectangle(0, 0, width, height),
                                                ImageLockMode.WriteOnly,
                                                PixelFormat.Format24bppRgb);

                unsafe
                {
                    byte* ptr = (byte*)bmpData.Scan0;
                    for (int i = 0; i < data.Length; i++)
                    {
                        Color color = HeatMapColor(data[i], min, max);
                        ptr[2] = color.R;
                        ptr[1] = color.G;
                        ptr[0] = color.B;
                        ptr += 3;
                    }
                }
                bmp.UnlockBits(bmpData);
                bmp.Save(path, ImageFormat.Png);
            }
        }

        private static Color HeatMapColor(byte value, byte min, byte max)
        {
            float ratio = (float)(value - min) / (max - min);
            return Color.FromArgb(
                (int)(255 * ratio),       // R
                (int)(255 * (1 - ratio)), // G
                0                         // B
            );
        }

        public static bool IsTileBlank(byte[][] tileData)
        {
            if (tileData == null || tileData.Length == 0) return true;

            int bandCount = tileData.Length;
            int pixelCount = tileData[0].Length;

            // 检查所有波段的所有像素是否一致
            byte baseValue = tileData[0][0];
            for (int b = 0; b < bandCount; b++)
            {
                byte[] data = tileData[b];
                if (data.Length != pixelCount) return false;

                foreach (byte val in data)
                {
                    if (val != baseValue) return false;
                }
            }
            return true;
        }

        public static List<GeoTileInfo> ReadGeoTileInfosFromFile(string geotiffSlitDir, bool isTheramlImg = true)
        {
            string tiffsplitsJsonFile = Path.Combine(geotiffSlitDir, isTheramlImg ? IR_GEO_TIFF_SPLIT_MAP_JSON_FILE_NAME : VS_GEO_TIFF_SPLIT_MAP_JSON_FILE_NAME);
            if (File.Exists(tiffsplitsJsonFile))
            {
                string content = File.ReadAllText(tiffsplitsJsonFile);
                if (!string.IsNullOrEmpty(content))
                {
                    return Newtonsoft.Json.JsonConvert.DeserializeObject<List<GeoTileInfo>>(content);
                }
            }
            return null;
        }

        public static GeoPoint GetGeoCoordinate(int pixelX, int pixelY,
            double[] geoTransform, string sourceProjection)
        {
            // 验证地理变换参数
            if (geoTransform == null || geoTransform.Length != 6)
                throw new ArgumentException("无效的地理变换参数");

            // 计算原始坐标
            double geoX = geoTransform[0] + pixelX * geoTransform[1] + pixelY * geoTransform[2];
            double geoY = geoTransform[3] + pixelX * geoTransform[4] + pixelY * geoTransform[5];

            // 创建坐标转换
            using (SpatialReference srcSrs = new SpatialReference(sourceProjection))
            using (SpatialReference destSrs = new SpatialReference(""))
            {
                destSrs.SetWellKnownGeogCS("WGS84");
                using (CoordinateTransformation ct = new CoordinateTransformation(srcSrs, destSrs))
                {
                    double[] output = new double[3];
                    ct.TransformPoint(output, geoX, geoY, 0);
                    return new GeoPoint
                    {
                        Lon = output[1],
                        Lat = output[0]
                    };
                }
            }
        }
    }

    public class FileGeoInfo
    {
        public Dataset Dataset { get; set; }
        public double[] Transform { get; set; }
        public CoordinateTransformation Ct { get; set; }
        public int XSize { get; set; }
        public int YSize { get; set; }
        public double MinX { get; set; }
        public double MaxX { get; set; }
        public double MinY { get; set; }
        public double MaxY { get; set; }
        public double PixelSizeX { get; set; }
        public double PixelSizeY { get; set; }
    }

    public class GridCoordinateSystem
    {
        public double OriginX { get; }
        public double OriginY { get; }
        public double PixelSizeX { get; }
        public double PixelSizeY { get; }

        public GridCoordinateSystem(double originX, double originY,
                                  double pixelSizeX, double pixelSizeY)
        {
            OriginX = originX;
            OriginY = originY;
            PixelSizeX = pixelSizeX;
            PixelSizeY = pixelSizeY;
        }
    }

}